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1.  Problem  of  Interest 

We  worked  on  a  fundamental  problem: 

•  Decompose  a  signal  into  a  small  set  of  decaying  complex  exponentials. 

Usually,  a  limited  sequence  {sjt}  of  noisy  samples  is  available  and  the  physics  of  the  problem  implies  that 
the  noiseless  exponentials  are  damped  sinusoids.  The  formal  signal  model  is 

Sk  =  '^arZ^~^,  (1) 

r=l 

where  the  coefficients  {ar}  and  knots  {zr}  are  2r  unknowns  to  be  determined.  This  problem  arises  in  a  wide 
range  of  disciplines,  including  nuclear  magnetic  resonance  [1],  speech  processing  [3],  and  system  identification 
[4].  Two  effective  solution  procedures  were  given  by  Prony  in  1795  [5]  and  S.Y.  Kung  et  aL  in  1983  [4]. 
Unfortunately,  Prony’s  method  is  relatively  unknown  axid  Kung’s  method  and  associated  proof  are  tightly 
bound  to  the  application  in  system  identification. 

1.1.  Our  Contributions 

We  developed  a  new  class  of  numerical  algorithms  for  solving  this  exponential  approximation  problem.  Our 
class  includes  both  methods  of  Prony  and  Kung  cis  special  cases.  We  gave  a  simple,  purely  linear  algebraic 
proof  on  why  our  new  approach  works. 

2.  Novel  Matrix  Pencil  Approach 


The  samples  {sk}  are  laid  out  as  entries  in  a  Hankel 

matrix: 

/5i 

52 

S3 

Sl 

52 

53 

S4 

•  5l+1 

jjKXL  ^ 

53 

54 

55 

•  5l-(-2 

Vsk 

5k+1 

5k+2  •  • 

*  Sk+L- 

where  K  >  R.  MatLab-like  notations  are  used  for  matrix  indexing,  e.g.,  the  colon  indicates  a  range  (hence 
^3^:  refers  to  the  third  row  of  A).  We  devised  a  new  matrix  pencil  procedure  for  computing  {zr},  after  which 
we  would  solve  a  special  (called  Yule- Walker)  set  of  linear  equations  to  find  {or}.  For  full  generality,  we 
introduce  two  arbitrary  but  nonsingular  transformations  F  and  G,  where  F  is  (k  -  1)  x  (k  —  1)  and  G  is 
L  X  L.  Define  two  new  (k  —  1)  x  l  matrices  and  by 

=  FHuk-i-G  and  =  FH2:k,.G. 

Consider  the  matrix  pencil  problem: 

=  CCt'ly. 

We  proved  that  solutions  exist  such  that  Y  is  a  L  x  R  matrix  of  linearly  independent  eigenvectors: 

cMy  =  (2) 

where  denotes  a  R  x  R  diagonal  matrix  of  eigenvalues  {C/}.  In  addition,  we  also  proved  that  {(;}  =  {zr}- 

Matrix  Pencil  Method. 

Step  1.  Choose  matrices  F  and  G. 

Step  2.  Solve  the  matrix  pencil:  =  C^^^YD^.  Set  {0}  = 

2.1.  Our  Contributions 

The  important  result  is  that  {0}  =  {zr},  so  that  eigenvalues  of  the  matrix  pencil  provide  the  required 
solution.  A  huge  benefit  is  the  freedom  to  choose  F  and  G.  Specific  choices  result  in  Prony ’s  and  Kung^s 
methods.  Other  choices  of  F  and  G  could  give  us  new  methods  with  other  desirable  characteristics. 
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3.  Full  Column  Rank  Case 

Both  Prony’s  and  Kung^s  methods  assume  full  column  rank,  i.e.,  L  =  R.  Hence  the  matrix  F  is  R  x  R  and 
nonsingular  and  (2)  becomes  Define  X  =  We  see  that  X  satisfies 

=  (3) 

The  matrix  equation  (3)  is  usually  overdetermined;  in  the  presence  of  noise,  X  may  be  solved  by  least-squares 
or  totcJ  least-squares  methods.  The  matrix  pencil  method  simplifies  to  an  eigenvalue  procedure. 

Eigenvalue  Method. 

Step  1.  Choose  matrices  F  and  G. 

Step  2.  Find  X  by  solving  the  linear  equations:  C^^^X  = 

Step  3.  Find  an  eigenvalue  decomposition  of  X:  X  =  FD^F“^.  Set  {0}  =  {^r}* 

3,1.  Prony’s  Procedure. 


Prony  further  assumed  that  K  =  R-h  1.  Hence  the  matrices 

and  are  R  x  R.  Choose  F  =  and 

G  =  So  =  /RxR  ^3^  simplifies  to  X  = 

Ct2l 

Note  that 

1 

0 

0 

0 

0 

1 
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0 

=  H2:R+l.l:RHn^l:R  = 
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0 
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0 
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0 

0 
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\-7o 

-7i 

_ frti 

-72 

•  •  •  -7r-2 

-7r-i  / 

where  the  {7,-}  satisfy  the  Yule- Walker  equation.  Hence  is  a  companion  matrix  to  the  polynomial  Pr{z), 
defined  by  Pk{z)  =  70  -f  7iz  H - h  7r_iz*^“^  +  z^, 

Prony’s  Method. 

Step  1.  Choose  F  =  and  G  = 

Step  2.  Solve  the  Yule- Walker  equation. 

Step  3.  Find  the  roots  {zr}  of  the  Prony  polynomial  Pk{z),  Set  {Q}  =  {zr}- 

3.2.  Kung’s  Procedure 

Start  by  computing  a  singular  value  decomposition  (SVD)  of  a  K  x  r  Hankel  matrix  H: 

Choose  F  =  /(K-i)x(K-i)  Q  ^  Then  and  =  £/2:k,:.  So 

Ui:K-i,:X  =  t/2:K,:-  An  important  advantage  of  Kung’s  method  over  Prony’s  procedure  is  that 
Pr{z)  is  not  explicitly  computed. 

Kung’s  Method. 

Step  1.  Choose  F  =  /(»<“1)x{k-i)  q  ^  yy;-!.  Compute  an  SVD  of  H\  H  = 

Step  2.  Solve  the  matrix  equation:  Uuk^i^-X  =  I72;k,:- 

Step  3.  Find  an  eigenvalue  decomposition  of  X:  X  =  FD^F“^.  Set  {Q}  =  {zr}* 

3.3.  Our  Contributions 

We  showed  that  our  matrix  pencil  scheme  includes  both  Prony’s  and  Kung’s  methods  as  special  cases.  So 
all  our  theoretical  results,  including  the  purely  linear  algebraic  proof  on  why  our  new  method  works,  cover 
these  two  procedures  too.  This  is  an  important  advance,  for  Kung’s  proof  [4]  can  be  difficult  to  digest. 


H  =  t/SF^. 
(3)  reduces  to 
the  polynomial 
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4.  New  Methods 

We  may  replace  the  expensive  SVD  by  other  techniques.  The  only  restriction  is  that  we  need  to  operate 
on  H  from  the  left  side.  An  example  of  a  much  cheaper,  and  yet  almost  cis  accurate,  technique  is  the  QR 
decomposition  (QRD).  Compute  a  QR  decomposition  of  a  k  x  r  Hankel  matrix  H:  H  =  QR,  Choose 
F  =  and  G  =  R-K  We  find  that  Gl‘l  =  gi:K-i,:  and  G^^l  =  Q2:k.:.  So  (3)  reduces  to 

QlrK-l,:-^  =  Q2:K,:  •  (4) 


Hankel  QBX)  Method. 

Step  1.  Choose  F  =  /(k-i)x(k-i)  and  G  = 

Step  2.  Solve  the  linear  equations:  Qv.k-i,:^  =  Q2:k,:- 

Step  3.  Find  an  eigenvalue  decomposition  of  X:  X  =  YD^Y~^,  Set  {0}  =  {zr}^ 

4.1.  Our  Contributions 

Our  new  Hankel  QRD  method  is  about  ten  times  faster  than  Kung’s  scheme  (also  known  as  the  Hankel  SVD 
method).  An  attractive  property  of  the  QRD  approach  is  that  it  is  easily  updatable  to  accommodate  new 
data,  which  is  not  so  for  an  SVD  technique. 

5.  Denoising 

Our  problem  is  particularly  sensitive  to  errors  in  the  samples  and  a  preliminary  noise  removal  procedure 
is  all-important.  We  assume  an  additive  error  model: 

If  =  i/W  4- ifW, 

where  and  denote  the  signal  and  noise  components  of  H,  respectively.  The  SVD  of  the  signal 
matrix  H  provides  a  convenient  and  numerically  stable  way  to  determine  the  noise  subspace.  Indeed,  it 
often  suffices  to  find  a  “gap”  in  the  entries  on  the  diagonal  of  E.  Unfortunately,  the  Hankel  structure  of  H 
is  almost  certainly  lost  in  Cadzow  [2]  proposed  to  iterate  signal  extraction  with  the  SVD  followed  by  a 
restoration  of  the  Hankel  structure  via  averaging  along  the  antidiagonals.  He  showed  that  under  relatively 
mild  assumptions  this  procedure  will  converge  and  reported  very  good  noise  removal  performance. 

5.1.  Our  Contributions 

First,  we  constructed  a  counter-example  for  Cadzow^s  method,  and  showed  that  it  will  converge  to  a  Hankel 
matrix  of  the  wrong  rank. 

Second,  we  proposed  a  new  approach  based  on  the  following  decomposition: 

H  =  (5) 

where  and  are  respectively  K  x  m  and  L  x  m  Vandemonde  matrices,  and  the  middle  matrix 

is  diagonal.  If  K  >  M,  L  >  M  and  all  the  {zr}  values  are  distinct,  then  the  Hankel  matrix  H  has  exact 
rank  M.  We  examined  algorithms  to  determine  the  factorization  (5)  for  a  given  Hankel  matrix.  It  seems 
reasonable  to  attempt  to  synthesize  another  Hankel  matrix  H  of  rank  R,  where  R  <  M,  from  zeroing  out 
(m  —  r)  of  the  {flr}  values.  However,  selecting  which  Or  to  set  to  zero  is  a  nontrivial  task.  Unlike  singular 
vectors,  the  columns  of  a  Vandemonde  matrix  are  not  orthonormal.  Hence,  it  is  not  always  appropriate 
to  select  the  R  largest  (in  magnitude)  Or’s  to  synthesize  a  good  approximation  to  H,  Our  work  is  still  in 
progress. 
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6.  Combined  Technique 

In  our  overall  computation,  we  chose  to  grossly  “overestimate”  the  rank  R  of  the  available  samples.  Just  as 
in  Prony’s  method,  we  would  pick  R  so  that  the  matrices  were  square.  Figure  1  shows  a  cloud  of  zeros  from 
the  resulting  polynomial  for  an  example  where  the  numerical  rank  of  the  signal  was  five,  but  a  65  x  65  noisy 
system  had  been  solved.  This  example  is  taken  from  Van  Huffel  [6];  the  signal-to-noise  ratio  is  approximately 
2.4  (not  dB). 
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Figure  1.  Prony^s  method  with  65  degrees  of  freedom  on  data  with  numerical  rank  5.  The  ‘-f* 
are  the  computed  z  values  while  the  exact  values  are  depicted  with 

With  the  common  cissumption  that  the  noise  space  is  approximately  orthogonal  to  the  signal  subspace,  we 
would  expect  that  a  few  components  will  still  fit  the  embedded  signal  while  most  other  zeros  serve  as  outlets 
to  fit  the  noise.  Our  numericcil  experiments  seem  to  confirm  this  conjecture.  Given  our  stability  criterion, 
it  is  logical  to  discard  those  zeros  (of  the  polynomial)  that  are  greater  than  one  in  magnitude  and  hence 
cannot  correspond  to  a  stable  component  of  the  signal.  Figure  2  shows  the  result. 


Figure  2.  left:  The  discrete  Fourier  spectrum  of  the  given  noisy  signal,  right:  The  spectrum 
after  removing  the  unstable  z  values  and  computing  the  new  samples  s.  The  dashed  curve  shows 
the  noiseless  spectrum  in  both  plots. 

6.1.  Our  Contributions 

There  is  a  significant  difference  between  our  approach  and  Prony’s  method:  we  fix  N,  the  number  of  samples, 
and  choose  R  =  [n/2J,  whereas  Prony  fixed  R,  the  rank  of  the  matrix,  and  chose  N  =  2r.  With  noise,  the 
Hankel  matrix  we  get  is  always  nonsingular.  Our  scheme  has  the  important  advantage  that  we  do  not  need 
to  decide  on  R  at  the  very  beginning.  Another  benefit  is  that  the  resultant  Yule- Walker  is  square  and  can 
be  solved  using  any  of  the  well  known  techniques. 
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